A  second  kind  integral  equation  formulation  is  presented  for  the  Dirichlet  problem  for  the 
Laplace  equation  in  two  dimensions,  with  the  boundary  conditions  specified  on  a  collection 
of  open  curves.  The  performance  of  the  obtained  apparatus  is  illustrated  with  several  nu¬ 
merical  examples.  The  formulation  is  a  simplification  of  the  equation  previously  constructed 
by  the  authors. 
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1  Introduction 


Integral  equations  have  been  one  of  principal  tools  for  the  numerical  solution  of  scattering  prob¬ 
lems  for  more  than  30  years,  both  in  the  Helmholtz  and  Maxwell  environments.  Historically, 
most  of  the  equations  used  have  been  of  the  first  kind,  since  numerical  instabilities  associated 
with  such  equations  have  not  been  critically  important  for  the  relatively  small-scale  problems 
that  could  be  handled  at  the  time. 

The  combination  of  improved  hardware  with  the  recent  progress  in  the  design  of  “fast” 
algorithms  has  changed  the  situation  dramatically.  Condition  numbers  of  systems  of  linear 
algebraic  equations  resulting  from  the  discretization  of  integral  equations  of  potential  theory 
have  become  critical,  and  the  simplest  way  to  limit  such  condition  numbers  is  by  starting  with 
second  kind  integral  equations.  Hence,  the  interest  in  reducing  scattering  problems  to  systems 
of  second  kind  integral  equations  on  the  boundaries  of  the  scatterers  has  been  rapidly  growing. 

During  the  last  several  years,  satisfactory  integral  equation  formulations  have  been  con¬ 
structed  in  both  acoustic  (Helmholtz  equation)  and  electromagnetic  (Maxwell’s  equations) 
environments,  whenever  all  of  the  scattering  surfaces  are  “closed”  (i.e.  scatterers  have  well- 
defined  interiors,  and  have  no  infinitely  thin  parts).  In  this  paper,  we  describe  a  second  kind 
integral  equation  formulation  for  the  Dirichlet  problem  for  the  Laplace  equation  with  boundary 
data  specified  on  a  collection  of  “open”  curves.  We  start  with  constructing  the  right  inverse  of 
the  single  layer  potential  operator  on  a  line  segment  via  simple  analytic  means;  then  we  apply 
such  operator  as  a  preconditioner  for  the  single  layer  potential  operator  on  the  curve  considered 
to  obtain  a  second  kind  integral  operator. 

Remark  1.1  In  a  recent  paper  [7],  the  authors  construct  a  somewhat  different  procedure  for 
the  solution  of  problems  of  the  classical  potential  theory  with  data  specified  on  a  collection  of 
open  surfaces.  While  the  approach  of  the  present  paper  is  very  similar  to  that  of  [7],  in  [7],  the 
single  layer  potential  is  used  to  precondition  the  quadruple  layer  potential  from  the  right;  here, 
the  quadruple  layer  potential  is  used  to  precondition  the  single  layer  potential  from  the  right. 
For  technical  reasons,  the  latter  leads  to  a  drastically  simplified  numerical  procedure  (and  also, 
requires  simpler  analysis);  hence,  this  sequel  to  [7]. 
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The  layout  of  the  paper  is  as  follows.  Section  2  contains  an  informal  description  of  the  pro¬ 
cedure.  In  Section  3,  the  necessary  mathematical  and  numerical  preliminaries  are  introduced. 
In  Sections  4,  we  present  the  principal  analytic  result  of  the  paper.  In  Section  5,  we  describe 
a  simple  numerical  implementation  of  the  scheme.  The  performance  of  the  algorithm  is  illus¬ 
trated  in  Section  6  with  several  numerical  examples.  Finally,  in  Section  7  we  discuss  several 
generalizations  of  the  approach. 


2  Informal  Description  of  the  Procedure 


In  this  section,  we  present  an  informal  description  of  the  procedure.  We  assume  that  7  : 
[— 1, 1]  -4  R2  is  a  sufficiently  smooth  “open”  (i.e.,  7(— 1)  ^  7(1))  curve  with  the  parametrization 


7(*)  =  7  (|  •(*+!)),  (1) 

where  L  is  the  total  arc  length  of  the  curve,  and  7  :  [0,  L]  -4  R2  is  the  same  curve  parametrized 
by  its  arc  length.  The  image  of  7  will  be  denoted  by  T.  We  consider  the  Dirichlet  problem  for 
the  Laplace  equation  in  two  dimensions,  with  the  boundary  conditions  specified  on  F,  i.e., 


Au  =  0  in  R2\r 
u  =  f  on  r. 


(2) 


This  problem  has  a  unique  bounded  solution  if  the  Dirichlet  data  /  is  sufficiently  smooth  (see, 
for  example,  [9]).  The  purpose  of  this  paper  is  to  reduce  the  problem  (2)  to  a  second  kind 
integral  equation  on  T. 

The  tools  of  the  classical  potential  theory  by  themselves  do  not  lead  to  such  an  integral 
equation.  Indeed,  the  standard  prescription  (see,  for  example,  [9])  is  to  represent  the  solution 
of  a  Dirichlet  problem  by  a  double  layer  potential,  and  the  solution  of  the  Neumann  problem 
by  a  single  layer  potential.  In  either  case,  the  behavior  of  the  singularity  near  the  boundary  is 
such  that  an  integral  equation  of  the  second  kind  on  T  is  obtained. 

However,  the  classical  procedure  critically  depends  on  T  being  a  closed  curve.  Indeed,  the 
potential  of  a  double  layer  on  the  curve  F  experiences  a  jump  when  T  is  crossed;  the  magnitude 


2 


of  the  jump  is  equal  to  the  density  of  the  double  layer  at  the  crossing  point.  This  poses  no 
problem  when  the  curve  is  a  closed  one,  since  the  potential  is  to  be  represented  on  only  one 
(inner  or  outer)  side  of  the  curve.  For  an  open  curve,  the  potential  has  to  be  represented  on 
both  sides  of  the  curve;  and  in  most  cases,  the  right-hand  side  /  (viewed  as  the  limiting  value  of 
the  solution  from  both  sides)  has  no  jump  across  F.  Thus,  an  attempt  to  represent  the  solution 
of  (2)  via  a  double  layer  potential  results  in  a  dipole  density  that  is  identically  equal  to  zero. 

One  could  attempt  to  represent  the  solution  of  (2)  by  a  charge  distribution  on  T.  The 
resulting  potential  is  continuous  across  T,  and  algorithms  of  this  type  have  been  constructed 
and  used  numerically  (see,  for  example,  [6]).  However,  the  resulting  integral  equation  is  of 
the  first  kind  (though,  fortunately,  with  a  logarithmically  singular  kernel),  with  all  the  usual 
numerical  disadvantages.  Another  option  is  to  use  the  quadruple  layer  potential  of  the  form 

R{cj)(x)  =  J  ^  a  (log \\x  -  7(011)  ’  °(t)dt,  (3) 

with  N(t)  the  unit  normal  to  T  at  7  (f);  the  resulting  equation  is  not  an  integral  equation  at 
all,  containing  a  part  that  is  actually  a  distribution.  In  engineering  literature,  such  objects  are 
known  as  “hypersingular  integral  equation” .  Satisfactory  procedures  have  been  constructed  for 
their  numerical  solution  (see,  for  example,  [3],  [10], [11]);  however,  these  are  not  as  simple  or  as 
stable  as  the  many  methods  available  for  the  solution  of  second  kind  integral  equations. 

This  paper  is  based  on  the  observation  that  when  the  curve  is  the  line  segment  I  =  [—1,1], 
the  right  inverse  of  the  single  layer  potential  operator  (denoted  by  Sj1)  can  be  constructed  by 
simple  analytic  means,  where  the  single  layer  potential  operator  Si  :  Ll[— 1,1]  — >  C[— 1,1]  is 
defined  by  the  formula 

Si{a)(x)  =  J  log  |a;  -  t\  ■  a(t)dt.  (4) 

Furthermore,  if  SJ1  is  used  as  a  preconditioner  for  the  single  layer  potential  operator  57  : 
Ll[— 1,1]  -4  C^M2)  on  T  defined  by  the  formula 

57(ct)(2;)  =  /^log|z-7(t)|  -a(t)dt,  (5) 
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i.e.,  the  solution  of  the  problem  (2)  is  represented  in  the  form 

u(x)  =  SyO  (6) 

then  the  resulting  boundary  integral  equation  is  of  the  second  kind. 

Remark  2.1  A  stable  second  kind  integral  equation  formulation  has  also  been  developed  for 
the  problem  (2)  in  [7].  Two  key  observations  used  in  [7]  are:  first,  the  product  of  the  quadruple 
layer  potential  operator  and  the  single  layer  potential  operator  is  a  second  kind  integral  operator 
for  the  case  of  a  closed  curve;  second,  the  case  of  a  line  segment  can  be  solved  analytically.  The 
integral  representation  for  the  solution  of  the  problem  (2)  in  [7]  is  of  the  form 

u(x)  =  Q1  o  5/  o  ( Q i  o  Si)~1(rj)(x),  (7) 

where  Qy  is  the  sum  of  a  quadruple  layer  potential  and  a  weighted  double  layer  potential  with 
the  weight  equal  to  the  curvature,  Si  is  the  single  layer  potential  operator  for  the  line  segment 
I  =  [—1, 1],  and  (Qi  o  5/)_1  is  (in  the  appropriate  sense)  the  right  inverse  of  Qi  o  5/.  The 
approach  of  this  paper  differs  from  that  of  [7]  in  that  the  roles  of  Q  and  S  are  interchanged, 
leading  to  a  simpler  scheme.  Indeed,  straightforward  analysis  shows  that  the  representation 
(6)  is  equivalent  to 

u(x)  =  S1oQJo(SioQI)~1(i1)(x).  (8) 

In  other  words,  the  solution  of  (2)  is  represented  by  a  single  layer  potential  on  F  preconditioned 
by  the  quadruple  layer  potential  for  the  line  segment  I,  with  a  further  preconditioning  by  the 
right  inverse  of  Si  o  Qi  to  eliminate  the  singularities  at  the  end  points. 

3  Analytical  Preliminaries 

In  this  section,  we  summarize  several  results  from  classical  and  numerical  analysis  to  be  used 
in  the  remainder  of  the  paper.  Detailed  references  are  given  in  the  text. 
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3.1  Chebyshev  Polynomials  and  Chebyshev  Approximation 

Chebyshev  polynomials  are  frequently  encountered  in  numerical  analysis.  As  is  well  known, 
Chebyshev  polynomials  of  the  first  kind  Tn  :  [-1, 1]  -*•  R  (n  >  0)  are  defined  by  the  formula 


Tn(x )  =  cos(narccos(x)), 


(9) 


and  are  orthogonal  with  respect  to  the  inner  product 
if, 9)  =  J :  fix)  •  g(x)  •  jl—dx. 


(10) 


Chebyshev  polynomials  of  the  second  kind  Un  :  [—1, 1]  -4  1  (n  >  0)  are  defined  by  the  formula 

sin  ((n  +  1)  arccos(x)) 


Un(x)- 

sm(arccos(x)J 

and  are  orthogonal  with  respect  to  the  inner  product 


(/,  g)  =  J  2  fix)  ■  gix)  ■  y/ 1  -  x2dx. 


The  Chebyshev  nodes  Xi  of  degree  N  are  the  zeros  of  T/v  defined  by  the  formula 

(2  i  +  1)7T 


Xi  =  cos 


2  N 


i  =  0, 1, . . . ,  N  -  1. 


(11) 


(12) 


(13) 


For  a  sufficiently  smooth  function  /  :  [ — 1, 1]  — >■  M,  its  Chebyshev  expansion  is  defined  by  the 
formula 

OO 

f(x)  =  Y^Ck-Tk(x),  (14) 

k= 0 

with  the  coefficients  Ck  given  by  the  formulae 

Co  =  ^  J1  fix)  •  T0{x)  •  (1  -  x2)~ldx,  (15) 

and 

ck  =  ^J  fix)  -Tkix)  ■  (1  -x2)~~%dx,  (16) 

for  all  A;  >  1.  We  will  also  denote  by  the  order  IV  —  1  Chebyshev  approximation  to 
the  function  /  on  the  interval  [—1,1],  i.e.,  the  (unique)  polynomial  of  order  IV  —  1  such  that 
Pfixi)  =  fixi)  for  all  i  =  0, 1, . . . ,  JV  —  1,  with  Xi  the  Chebyshev  nodes  defined  by  (13). 

The  following  lemma  provides  an  error  estimate  for  the  Chebyshev  approximation  (see,  for 
example,  [4]). 
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Lemma  3.1  If  f  £  Ck[— 1,1]  (i.e.,  f  has  k  continuous  derivatives  on  the  interval  [ —  1 , 1}), 
then  for  any  x  £  [— 1, 1], 

|Pf(x)  -  /(*)|  =  O  .  (17) 

In  particular,  if  f  is  infinitely  differentiable,  then  the  Chebyshev  approximation  converges  su- 
peralgebraically  (i.e.,  faster  than  any  finite  power  of  l/N  as  N  —>•  ooj. 

3.2  Miscellaneous  Results 


In  this  section,  we  collect  several  results  from  classical  analysis  to  be  used  subsequently.  Lemma 
3.2  lists  two  standard  definite  integrals;  both  can  be  found  (in  a  somewhat  different  form)  in 
[5].  Lemma  3.3  states  a  standard  fact  from  classical  potential  theory;  it  can  be  found  in  [9]. 
Finally,  Lemma  3.4  states  that  if  the  curve  7  is  sufficiently  smooth,  then  the  restriction  of  the 
kernel  of  the  operator  57  —  5/  on  T  is  also  smooth  (see  (4),  (5)  for  the  definitions  of  Si  and 
S7). 


Lemma  3.2  For  any  x  £  (—1, 1), 


and 


r1  1 

J  1  log  \x-t\-  -y==dt  =  -7T  •  log  2, 


rl  1  _ 

p.v.  /  - -  -Un-i{t)  •  VI  -t2dt  =  71  ■  Tn(x), 

J —1  x  t 


(18) 


(19) 


for  any  n  >  1. 


Lemma  3.3  Suppose  that  7  :  [—  1, 1]  — >  M2  is  a  sufficiently  smooth  open  regular  curve  with  the 
parametrization  (1),  and  that  the  function  a  £  Ll[— 1,1]  satisfies  the  condition 

r»l 

a(t)dt  =  0.  (20) 


/: 

Then  the  function  u  :  R2  — >•  R  defined  by  the  formula 
u(x)  =  J  log  \x  —  7(i)|  •  a(t)dt 


(21) 


is  bounded  in  M2 . 
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Lemma  3.4  Suppose  that  7  G  Ck+l  [0,  L\  (k  >  1)  is  an  open  regular  curve  parametrized  by 
its  arc  length  in  R2 .  Suppose  further  that  the  function  r  :  [0,  L]  x  [0,  L]  — >  M  is  defined  by  the 
formula 


r(x,  t)  =  < 


f  log  |7(x)  -  7(i)|  -  log  \x-t\,  x^t 
(0,  x  =  t. 

Then  r  G  Ck([0,L]  x  [0,L]). 


(22) 


Proof.  Since  7  is  parametrized  by  its  arc  length,  we  have 
W(x)\  =  1, 

for  all  x  G  [0,  L).  Combining  (22),  (23),  we  observe  that 
r(x,t)  —  log  \h(x,t)\, 

where  the  function  h  :  [0,  L]  x  [0,  L]  ->  M2  is  defined  by  the  formula 

7(z)  -  7(f) 


h(x,  t)  =  { 


x  —  f 


-,  x  ^t 


(23) 


(24) 


(25) 


l  7' (a;),  x  =  t. 

Obviously,  h  is  k  times  continuously  differentiable  for  7  €  Ck+1[0,L\  by  Taylor’s  Theorem. 
Furthermore,  since  7(2)  ^  7(f)  if  x  ^  t,  and  |7/(m)|  =  1  for  all  x  G  [0,  L\,  we  have 

\h(x,  t)|  7^  0  for  all  (x,t)  G  [0,L]  x  [0,L].  (26) 

Therefore,  the  function  r  =  log \h\  is  also  k  times  continuously  differentiable  in  [0,  L]  x  [0,  L). 


□ 


4  Analytical  Apparatus 

4.1  Right  Inverse  of  the  Single  Layer  Potential  Operator  on  the  Line  Seg¬ 
ment 

The  purpose  of  this  section  is  Theorem  4.2,  providing  the  right  inverse  of  the  single  layer 
potential  operator  on  the  line  segment  I  —  [—1, 1].  The  construction  is  based  on  an  elementary 
integral  identity  stated  in  Lemma  4.1. 
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Lemma  4.1  For  any  x  G  (—1, 1), 

J  ^  log  |x  -  t\  ■  J^pdt  =  -7r  ■  log  2  •  T0(x),  (27) 

and 

[  \og\x-t\--^=^dt  =  -^-Tn{x),  (28) 

J-i  Vl  -t2  n 

for  any  n  >  1. 

Proof.  (27)  directly  follows  from  the  combination  of  the  identity  (18)  and  the  fact  that 
Tq(x)  =  1  for  all  x  G  [—1, 1].  To  prove  (28),  we  integrate  by  parts  once,  obtaining 

f  log  |x  —  f|  •  -^=L=dt  =  — — p.v.  f  -}—-Un-i(t)-Vl-t2dt.  (29) 

J- 1  Vl  -  t2  n  x  -  t 

Now,  (28)  follows  from  the  combination  of  (29),  (19).  □ 


Theorem  4.2  Suppose  that  the  linear  operator  S  :  C[—  1, 1]  -»  L1[— 1, 1]  is  defined  by  its  action 
on  the  functions  Tn  (n  >  0)  via  the  formula 

1  T0(x) 

=  ,  n  =  u 

(30) 


S(Tn)(x)  = 


7T-log2  v/l  -  X2’ 

n  Tn(x) 


n  =  0 
n  >  0. 


7T  y/1  -  x2  ’ 

Suppose  further  that  the  operator  Si  :  L1[— 1, 1]  -»  C[—  1, 1]  is  defined  by  the  formula 


Si(a)(x)  =  J  ^ 


log  \x  -t\-  a(t)dt. 


(31) 


Then 


Si  o  5  =  T,  (32) 

with  T  the  identity  operator.  In  other  words,  S  is  the  right  inverse  of  Si  on  the  space  of 
continuous  functions. 
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Proof.  Since  Tn  (n  >  0)  form  a  basis  for  the  space  C[—  1, 1],  and  the  operators  Si,  S  are  linear, 
we  only  need  to  prove  that  the  identity 


Si  o  S(Tn)(x)  =  Tn(x) 


(33) 


holds  for  all  n  >  0.  Substituting  (30)  into  (31)  we  obtain 


Si  o  S(Tn)(x)  = 


-  ■  [  log  \x  - 

7T 


t\  ■  -~^==dt, 

1  VT^ 


j  log  |x  t\ 


Tn(t) 


Vl-^2 

which  directly  follows  from  Lemma  4.1 


dt  — 


n 


n  =  0 
n  >  0. 


Combining  (33),  (34),  we  observe  that  it  suffices  to  prove  the  identity 

7r  •  log  2  •  !Zo(x),  n  —  0 
n  >  0, 


(34) 


(35) 


□ 


4.2  Second  Kind  Integral  Equation  Formulation 

In  this  section,  we  reduce  Problem  (2)  to  an  integral  equation  of  the  second  kind  on  the 
curve  T;  the  results  are  summarized  in  Theorem  4.4.  We  start  with  defining  the  operator 
S7  :  C[—  1, 1]  — >  C(M2)  via  the  formula 


57(o-)(z)  =  S1oS(a)(z), 


(36) 


with  S7,  S  defined  by  (5),  (30),  respectively.  Combining  (36)  with  Theorem  4.2,  we  easily  see 
that  for  arbitrary  smooth  a  :  [—1, 1]  ->  E  and  'y(x)  €  T, 

57(<r)( 7(2))  =  5/o  S(a){x)  +  (57  -  5/)  o  5(ct)(7(®)) 

=  cr(a>)  +  (57  -  Si)  o  S(a)( j(x)), 

and  the  following  theorem  shows  that  the  operator  P7  =  (S7  —  Si)  o  5  is  compact. 


(37) 
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Theorem  4.3  Suppose  that  7  :  [—1, 1]  -4  R2  is  a  sufficiently  smooth  open  regular  curve  with 
the  parametrization  (1).  Suppose  further  that  the  operator  P 7  :  C[—l,  1]  — »  C[— 1,1]  is  defined 
by  the  formula 

P y(a)(x)  =  (S1  -  Sj)  o  S{a)( 7(0;)) 

=  y  ^(l°g  |7(X)  —  7(01  —  log  I1  ~  *1)  ’  S(a)(t)dt,  ^ 

with  Sy,  Si,  S  defined  by  (5),  (31),  (30),  respectively.  Then  P7  is  compact. 

Proof.  By  Lemma  3.4,  the  function  r  :  [-1, 1]  x  [—1, 1]  ->  M  defined  by  the  formula 


r(x,  t)  =  log  |7(:e)  -  7(t)|  -  log  \x  -  t\ 


(39) 


is  k  times  continuously  differentiable  for  any  7  €  Cfc+1[— 1, 1].  Obviously,  if  r  is  expanded  into 
a  double  Chebyshev  series 

OO  OO 

?(*.<)  =  ££  P-mnTm  (x)Tn(t),  (40) 

m=0n=0 

then  there  exists  a  positive  number  C  such  that 


\K„ 


< 


mk  ■  nk 


for  any  m  >  0,  n  >  0.  Now,  for  any  N  >  0,  we  will  define  the  operator  P/v  :  C[—  1, 1] 
by  the  formula 


(41) 

C[-l,l] 


PN(cr){x)  =  j  rN(x,t)  ■  S(o)dt, 

with  the  function  rjv  :  [-1. 1]  x  [— 1, 1]  ->  R  defined  by  the  formula 

N  N 


(42) 


rN(x,t)  =  YJY,Kmnrm(x)Tn{t).  (43) 

m= 0  n=0 

Obviously,  Pn  is  a  compact  operator  since  its  range  is  of  finite  dimensionality.  Furthermore, 
PN  converges  to  P7  as  TV  ->  oo  by  (41).  Hence,  P7  is  also  a  compact  operator.  □ 
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We  will  represent  the  solution  of  Problem  (2)  via  the  formula 


u(x )  =  Sy(a)(x)  +  A 

=  [  log|x-7(t)|  -S(cr)(t)dt  +  A, 


(44) 


where  A  is  a  real  constant  to  be  determined.  Combining  Lemma  3.3  and  Theorem  4.3,  we 
obtain  the  principal  result  of  this  paper. 


Theorem  4.4  Suppose  that  7  :  [—1, 1]  R2  is  a  sufficiently  smooth  open  regular  curve  with 
the  parametrization  (1),  and  that  the  function  f  :  [—1, 1]  — ►  M  is  continuously  differentiable. 
Suppose  further  that  the  continuous  function  a  :  [—  1, 1]  — >  R  and  the  coefficient  A  satisfy  the 
equations 


a{x)  +  Py(cr)(x)  =  f{x)  -  A, 

/'^•Trb^0' 


(45) 

(46) 


with  Py  defined  in  (38) .  Then  the  function  u  :  R2  -»  M  defined  by  (44)  is  bounded  in  M2  and  is 
the  solution  of  the  problem 


A u  =  0  in  K2\r 

<  (47) 

u  =  f  on  T. 


Remark  4.1  Obviously,  the  purpose  of  the  constant  A  in  the  above  theorem  is  to  ensure  the 
boundedness  of  the  solution  u  of  (2).  In  certain  physical  situations,  the  potentials  of  interest 
are  not  bounded  at  infinity,  but  rather  grow  logarithmically.  In  such  cases,  the  solution  to  (2) 
assumes  the  form 


u(x)  =  Sy(cr)(x), 

with  a  satisfying  the  integral  equation 


(48) 


<r(x)  +  Py(a)(x)  =  f(x). 


(49) 
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5  Numerical  Algorithm 


In  this  section,  we  construct  a  rudimentary  numerical  algorithm  for  the  solution  of  the  Dirichlet 
problem  (47)  via  the  equations  (45)  -  (46).  Since  the  construction  of  the  matrix  and  the  solver 
of  the  resulting  linear  system  are  direct,  the  algorithm  requires  0(N3)  work  and  0(N2)  storage, 
with  N  the  number  of  nodes  on  the  boundary.  While  standard  acceleration  techniques  (such  as 
the  Fast  Multipole  Method,  etc.)  could  be  used  to  improve  these  estimates,  no  such  acceleration 
was  performed,  since  the  purpose  of  this  section  (as  well  as  the  following  one)  is  to  demonstrate 
the  stability  of  the  integral  formulation  and  the  convergence  rate  of  a  very  simple  discretization 
scheme. 

By  Theorem  4.4,  the  equations  to  be  solved  axe  (45)  -  (46),  where  the  unknowns  are  the 
function  a  and  the  real  number  A.  To  solve  (45)  -  (46)  numerically,  we  discretize  the  boundary 
into  N  Chebyshev  nodes  and  approximate  the  unknown  density  a  by  a  finite  Chebyshev  series 
of  the  first  kind, 


N-l 


cr(t)~^Ck-Tk(t), 

fc= 0 


(50) 


with  the  coefficients  Ck  (fc  =  0, . . . ,  N  —  1)  to  be  determined.  In  order  to  discretize  (45),  we 
start  with  observing  that  by  (29),  the  action  of  the  operator  S  on  the  function  a  is  described 
via  the  formula 


S(ct)(x)  = 


N- 1 


Tk(x), 
k= o 


(51) 


where  the  coefficients  Bk  (A;  =  0, . . . ,  TV  -  1)  are  given  by  the  formulae 
B0  = - :  X 

< 


7T  log  2  ’ 


(52) 


Bk  =  —  1  <  k  <  N  -  1. 

7T 


Next,  we  approximate  the  kernel  r(x,  t)  (see  (40))  of  the  operator  S7  —  Si  with  an  expression 
of  the  form 

N-1N-1 

r(x,t )  ~  53  Kij  ■  Ti(x)  ■  Tj(t).  (53) 

t=0  j= o 
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Cleaxly,  the  coefficients  Kij  have  to  be  determined  numerically,  since  the  curve  F  is  user- 
specified,  and  is  unlikely  to  have  a  convenient  analytical  expression.  Thus,  we  obtain  the 
coefficients  Kij  by  first  constructing  the  N  x  N  matrix  R  -  ( r(xi,tj ))  ( i,j  =  0, 1, . . . ,  N  -  1) 
with  Xi,tj  the  Chebyshev  nodes  defined  by  (13)  then  converting  R  into  the  matrix  K  =  (Kij) 
(i,  j  =  0, 1, . . . ,  N  —  1)  by  the  formula 

K  =  U-RUT,  (54) 

with  N  x  N  matrix  U  =  (Uij)  defined  by  the  formula 

u0j  =  ^7  •  To (xj),  j  =  0, 1, . . . ,  N  -  1, 

Ti(xj),  i=l,...,N-l,  j  =  0,1, . . .  ,N  —  1, 

Finally,  we  approximate  the  prescribed  Dirichlet  data  /  by  its  Chebyshev  approximation  of 


3  N 

Uij  =  — 

N 


(55) 


order  N  —  1 


JV-l 

fc= 0 


(56) 


where  the  coefficients  /*,  can  be  obtained  by  first  evaluating  /  at  Chebyshev  nodes  xu  then 
applying  to  it  the  matrix  U  defined  by  (55),  i.e., 

N-l 

Uki  •  f(Xi).  (57) 

t=0 

Combining  (51),  (53),  (56),  we  discretize  (45)  into  the  equation 


A  ■ 


(  c°  ) 

(  1  ^ 

f  \ 

Cl 

+  A- 

0 

= 

A 

\  Cn-i  J 

U  ) 

\  In-i  j 

(58) 


with  N  x  N  matrix  A  defined  by  the  formula 

A  =  IN  +  K  ■  B,  (59) 

with  7jv  the  N  x  N  identity  matrix,  and  B  the  diagonal  matrix  defined  by  the  formula 

Bij  =  Bi  •  Sij.  (60) 
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Furthermore,  (44)  leads  to  the  equation 


Co  =  0, 


(61) 


Finally,  (58)  and  (61)  together  form  a  linear  system  of  dimension  N  +  1  to  be  solved. 

Remark  5.1  The  generalization  of  the  above  scheme  to  the  case  of  several  disjoint  open  curves 
is  straightforward,  and  has  been  implemented  by  the  authors  (see  Example  4  in  Section  6). 

6  Numerical  Examples 

A  FORTRAN  code  has  been  written  implementing  the  algorithm  described  in  the  preceding 
section.  In  this  section,  we  demonstrate  the  performance  of  the  scheme  with  several  numerical 
examples.  We  consider  the  problem  in  electrostatics:  the  boundary  is  made  of  conductor  and 
grounded,  the  electric  field  incident  on  the  boundary  is  generated  by  the  sources  outside  the 
boundary.  For  these  examples,  we  plot  the  equipotential  lines  of  the  total  field  and  present 
tables  showing  the  convergence  rate  of  the  algorithm. 

Remark  6.1  In  the  examples  below,  the  problems  to  be  solved  via  the  procedure  of  the  preced¬ 
ing  section  have  no  simple  analytical  solution.  Thus,  we  tested  the  accuracy  of  our  procedure 
by  evaluating  our  solution  via  the  formula  (44)  at  a  large  number  M  of  nodes  on  the  boundary 
T  (in  our  experiments,  we  always  used  M  =  2000),  and  comparing  it  with  the  analytically 
evaluated  right-hand  side.  We  did  not  need  to  verify  the  fact  that  our  solutions  satisfy  the 
Laplace  equation,  since  this  follows  directly  from  the  representation  (44). 

In  each  of  those  tables,  the  first  column  contains  the  total  number  N  of  nodes  in  the 
discretization  of  each  curve.  The  second  column  contains  the  condition  number  of  the  linear 
system.  The  third  column  contains  the  relative  L2  error  of  the  numerical  solution  as  compared 
with  the  analytically  evaluated  Dirichlet  data  on  the  boundary.  The  fourth  column  contains 
the  maximum  absolute  error  on  the  boundary.  In  the  last  two  columns,  we  list  the  errors  of  the 
numerical  solution  as  compared  with  the  numerical  solution  with  twice  the  number  of  nodes, 
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Table  1:  Numerical  results  for  Example  1. 


N 

K 

E'\T) 

e°°(  r) 

E'\u) 

E°°{u) 

8 

0.200 E  +  01 

0.703E  -  01 

0.17822  +  00 

0.29622  -  02 

0.52822  -  02 

16 

0.222E  4-  01 

0.759E  -  02 

0.21222  -  01 

0.64122  -  04 

0.11425  —  03 

32 

0.212E  +  01 

0.165E-03 

0.48622  -  03 

0.55622  -  07 

0.99122  -  07 

64 

0.206E  +  01 

0.14712  -  06 

0.44622  -  06 

0.83522  -  13 

0.15022  —  12 

128 

0.203J5  +  01 

0.225E  -  12 

0.69022  -  12 

0.35522  -  15 

0.22222  -  14 

256 

0.20225  +  01 

0.93522  -  15 

0.21422  -  13 

0.34322  -  15 

0.20022  -  14 

where  the  solution  is  evaluated  at  1000  equispaced  points  on  a  circle  of  radius  3.3  centered  at 
the  origin;  the  fifth  column  contains  the  relative  L 2  error,  and  the  sixth  column  contains  the 
maximum  absolute  error. 

Example  1:  In  this  example,  the  boundary  is  the  line  segment  parametrized  by  the  formula 


f  x(t )  =  t 

1  y(t )  =  -0.2 


-  1  <  t  <  1. 


(62) 


The  Dirichlet  data  axe  generated  by  a  unit  charge  at  (0,0).  The  numerical  results  are  shown 
in  Table  1.  The  source,  curve  and  equipotential  lines  are  plotted  in  Figure  1. 


Figure  1:  Source,  curve,  and  equipotential  lines  for  Example  1. 
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Table  2:  Numerical  results  for  Example  2. 


N 

K 

E'\Y) 

E°°(r) 

E\u) 

E°°(u) 

32 

0.195E  +  01 

0.271  E  -  01 

0.8641?  -  01 

0.6581?  -  02 

0.4691?  -  02 

64 

0.187E  +  01 

0.2401?  -  02 

0.8471?  -  02 

0.1461?  —  03 

0.1041?  —  03 

128 

0.182E  +  01 

0.4221?  -  04 

0.1571?  —  03 

0.1351?  —  06 

0.9551?  -  07 

256 

0.1791?  +  01 

0.3071?  -  07 

0.1171?  —  06 

0.2451?  -  12 

0.1731?  —  12 

512 

0.1781?  +  01 

0.4311?  -  13 

0.1601?  —  12 

0.9711?  —  15 

0.133 E  -  14 

1024 

0.1771?  +  01 

0.3041?  -  14 

0.4501?  -  13 

0.941  E  -  15 

0.1221?  —  14 

Example  2:  In  this  example,  the  boundary  is  a  sinusoidal  arc  parametrized  by  the  formula 

(63) 


{ 


x(t) 

y(t) 


0.5 1 
cos  (t) 


37 r  37r 

—  <  t  <  — . 
2  ~  -  2 


The  Dirichlet  data  are  generated  by  one  positive  charge  of  unit  strength  at  (0, 1.5)  and  another 
negative  charge  of  unit  strength  at  (0,0).  The  numerical  results  are  shown  in  Table  2.  The 
sources,  curve  and  equipotential  lines  are  plotted  in  Figure  2. 


Figure  2:  Sources,  curve,  and  equipotential  lines  for  Example  2. 


Example  3:  In  this  example,  the  boundary  is  a  spiral  parametrized  by  the  formula 
x(t)  =  f  cos(3.37rt)  —  0.1 


y(t)  =  tsin(3.37rt) 


0.2  <  t  <  3.2. 


(64) 
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Table  3:  Numerical  results  for  Example  3. 


N 

K 

E2(T) 

e°°(  r) 

E2(u) 

E°°(u) 

32 

0.704 E  +  03 

0.594F  -  01 

0.125F  +  00 

0.23325  4-  00 

0.68525  -  01 

64 

0.65715  +  02 

0.108F  -  02 

0.665E  -  02 

0.41725  -  02 

0.201  E  -  02 

128 

0.523E  +  02 

0.904F  -  04 

0.653 E  -  03 

0.10125  -  03 

0.575 E  -  04 

256 

0.394E  +  02 

0.213F  -  05 

0.18315  —  04 

0.17925  -  06 

0.12525  —  06 

512 

0.279F?  +  02 

0.313F  -  08 

0.27215  -  07 

0.15625  -  11 

0.12325  —  11 

1024 

0.196F  +  02 

0. 18415  —  13 

0.14725  -  12 

0.21125  —  13 

0.93315  - 14 

The  Dirichlet  data  are  generated  by  a  unit  charge  at  (0,0).  The  numerical  results  are  shown 
in  Table  3.  The  source,  curve  and  equipotential  lines  are  plotted  in  Figure  3. 


Figure  3:  Source,  curve,  and  equipotential  lines  for  Example  3. 


Example  4:  In  this  example,  we  consider  the  case  of  several  open  curves.  The  boundary 
consists  of  three  elliptic  arcs  parametrized  by  the  formulae 


f  *i(<) 
l  Vi(t) 

f  *2  (t) 
I  V2  (t) 


=  — tcos(3.37rt)  —  1.45 
=  — tsin(3.37rt)  -f  0.55 


0.2  <  t  <  1.2, 


=  £cos(3.37rf)  —  0.1 
=  tsin(3.37rt)  —  1.2 


0.2  <t<  1.2, 


(65) 

(66) 
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